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Instability of spatial patterns and its ambiguous impact on species diversity 
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Self- arrangement of individuals into spatial patterns often accompanies and promotes species 
diversity in ecological systems. Here, we investigate pattern formation arising from cyclic dominance 
of three species, operating near a bifurcation point. In its vicinity, an Eckhaus instability occurs, 
leading to convectively unstable "blurred" patterns. At the bifurcation point, stochastic effects 
dominate and induce counterintuitive effects on diversity: Large patterns, emerging for medium 
values of individuals' mobility, lead to rapid species extinction, while small patterns (low mobility) 
promote diversity, and high mobilities render spatial structures irrelevant. We provide a quantitative 
analysis of these phenomena, employing a complex Ginzburg- Landau equation. 
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Biodiverse ecosystems comprise complex interactions 
of a large number of individuals and species leading to 
rich spatio-temporal community structures f]\. Much 
efforts in theoretical and biological physics are cur- 
rently devoted to qualitative and quantitative under- 
standing of basic mechanisms that maintain their di- 
versity. Hereby, within exemplary models, the forma- 
tion of dynamic spatial patterns has been identified as 
a key promoter [2 B S [E|- In particular, crucial influ- 
ence of self-organized patterns on biodiversity has been 
demonstrated in recent experimental studies p , employ- 
ing three bacterial strains that display cyclic competi- 
tion. The latter is metaphorically described by the game 
'rock-paper-scissors' where rock smashes scissors, scis- 
sors cut paper, and paper wraps rock in turn. Such 
non-hierarchichal dynamics has also been found in, e.g., 
lizard populations in California [7] and coral reef inverte- 
brates 0. For the three bacterial strains, and for low mi- 
crobes' motility, cyclic dominance leads to stable coexis- 
tence of all three strains through self-formation of spatial 
patterns [6]. In contrast, stirring the system, as can also 
result from high motilities of the individuals, destroys the 
spatial structures and the strain coexistence 0. 

Here, from theoretical studies, we show that cyclic 
competition of species can lead to highly non-trivial spa- 
tial patterns as well as counterintuitive effects on biodi- 
versity. Investigating cyclic dynamics near a (degenerate) 
bifurcation, we find that pattern formation is only weak 
and mainly influenced by stochastic effects. Namely, in 
a prototypical model where three species compete cycli- 
cally, and the deterministic rate equations (RE) for the 
densities' time evolution predict weakly instable coexis- 
tence (near or at neutral stability), we demonstrate the 
generic formation of convectively instable spiral waves. 
The instability of the latter appears as a "blurring" of the 
spirals and is the stronger the closer the RE operate near 
the bifurcation point, i.e. near the parameter point where 
the RE predict neutral stability. Consequently, patterns 



take shape only weakly and allow for major influence of 
stochasticity. This effect is most pronounced at the bi- 
furcation point; there, noise remains as the only source 
of pattern formation. Furthermore, at the bifurcation 
point, we uncover a counterintuitive destabilizing effect of 
patterns on the stability of coexistence. Similar to what 
has been found in other studies, see Refs. [1, 0, 0, [13, [HI 
and references therein, for low individuals' mobility, the 
size of the emerging spatial structures lies much below the 
size of the ecosystem and allows for stable coexistence of 
all species. However, we show that a mobility exceeding 
a first critical value leads to rapid extinction of species 
with only one surviving. This scenario is connected to 
weakly formed patterns that span over the whole sys- 
tem, revealing the antagonistic character of such large 
patterns on biodiversity. Only when mobility lies above 
a second threshold, a third, effectively well-mixed regime 
emerges. There, patterns and spatial correlations do not 
form and the fate of species diversity is determined by 
the character of the species competitions alone. 

We design a spatial and stochastic model of cyclically 
competing populations in the following way. Individu- 
als of three species A, and C populate a square lat- 
tice, such that each site is occupied by one individual or 
left empty. Each agent can interact with its four near- 
est neighbors by either exchanging positions at a rate e 
(individuals are mobile), or by cyclic competition. The 
latter interactions are described in the language of chem- 
ical reactions; generically, we consider: 

AB AA, AB ^ A0, A0 AA, 

BC BB, BC ^ B0, B0 BB, 

CA CC, (1) CA C0, (2) C0 CC. (3) 

Reactions ([T]) describe consumption of an individual by 
another of a predominant species, and immediate repro- 
duction of the latter. Cyclic dominance appears as A 
consumes B and reproduces, while B preys on C and 
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FIG. 1: Snapshots of the biodiverse state for D = 1 x 10~^. 

(a) , For large rates 7, entangled and stable spiral waves form. 

(b) , A convective (Eckhaus) instability occurs at 7^; ^ 2; be- 
low this value, the spiral patterns blur, (c), At the bifurcation 
point 7 = 0, only very weak spatial modulations emerge; we 
have amplified them by a factor two for better visibility. The 
snapshots stem from numerical solution of the SPDE (|5)) with 
initially homogeneous densities a — h — c— l/A. 



C feeds on A in turn. We consider the most symmetric 
version where all species are symmetric, and fix the time- 
scale by setting the rates for these reactions to one. Re- 
actions ([2]) encode solely consumption, leaving an empty 
site 0. These reactions occur at a rate 7, and are de- 
coupled from reproduction, Eqs. (|3|), which happens at 
a rate /i. Note that reactions ([1]) and (|3|) describe two 
different mechanisms of reproduction, both of which are 
important for ecological systems: In ([1]), an individual re- 
produces when having consumed a prey, due to thereby 
increased fitness. In contrast, in reactions ([3|) reproduc- 
tion depends solely on the availability of empty space. 

The RE encode the deterministic behavior of these re- 
actions in the well-mixed scenario. Dependent on the 
type of reactions, they can yield stable, unstable, or neu- 
trally stable diversity. For the above defined model, with 
X = {a^b^c) denoting the densities of the three distinct 
species and p = a -\- b -\- c the total density, the RE read, 

dtXi = Xi - p) + Xi^i - (1 + 7)^i+2] = Ai . (4) 

Hereby, the indices are understood as modulo 3. These 
equations have been analyzed by May and Leonard [l2[ 
and also appear in theoretical descriptions of the 
Kuppers-Lortz instability [l^, Generically, Eqs. (|4|) 
possess an unstable internal fixed point as well as hetero- 
clinic orbits, where the system cycles between states with 
nearly only one species, leading to rapid extinction of all 
but one species when stochastic effects are included. A 
degenerate bifurcation emerges at /i = 7 = 0, i.e. when 
only the reactions ^ are present. In this situation, neu- 
trally stable orbits surround an internal fixed point, being 
neutrally stable as well [H, [l5[ . Finite-size fluctuations 
invalidate the neutral stability and induce extinction of 
two species after a characteristic time T which is pro- 
portional to the system size 



15|. In the following. 



this bifurcation point; for illustration, we consider equal 
selection and reproduction rates 7 = /i. 

In the spatial model, mobility of individuals, stemming 
from random local exchanges, leads to their diffusion on 
the lattice with a diffusion constant D = 2eN~^ ll, 
Employing a continuum limit of large systems, N ^ 00 
with diffusivity D kept fixed, the stochastic spatial sys- 
tem becomes describable by stochastic partial differential 
equations (SPDE), see References [l^; they read 



dtXi = DV^Xi^Ai^'^Cij^j 



(5) 



In these equations, a term D\/^ describes individuals' 
diffusion on the lattice; it reveals that the size of possible 
spatial structures is proportional to 13]. The non- 



linear terms Ai coincide with those of the RE (|4|) and 
take interactions into account. Noise terms, stemming 
from a system-size expansion 18, [l9| are added which 
scale as the square root of N. In Eqs. ([5]), denote un- 
corr elated gaussian noise, and correlations are accounted 
for by CC^. The matrix C is thereby not uniquely deter- 
mined. While CC^ is symmetric, C appears, in general, 
asymmetric without physical significance; we choose 



C = 



2^7 
6N V2V7 



-V3V3 + 27 


73 V3 + 27 



-V3 + 27 
2V3 + 27 
-V3 + 27, 



(6) 



We have numerically solved the SPDE ([5]) using open 
software from the XmdS project [13]. Snapshots of the 
resulting steady states, for low diffusivity and rates 7 
approaching the bifurcation point 7 = 0, are shown in 
Fig. [H All three species coexist in a stable manner, 
through the formation of spatial patterns. Large values 
of 7 (see 7 = 5 in Fig. [1] (a)) yield a dominant contribu- 
tion of the selection and reproduction events ([2]),(|3|), and 
entangled rotating spirals form, similar to the structures 
reported in Refs. [5|, ll7[. At lower values, these spirals 
appear blurred, while their wavelength increases. The 
blurring intensifies upon lowering 7; at the bifurcation 
point (7 = 0) patterns take shape only extremely weakly 
and appear to be of predominantly stochastic nature. 

These observations can be analytically understood 
by employing a complex Ginzburg-Landau equation 
(CGLE). Although the RE g]) operate in a three- 
dimensional phase space, spanned by the densities a, 6, 
and c, trajectories quickly relax to a two-dimensional in- 
variant manifold [l7[. On the latter, expanding the RE 
to third order around the unstable reactive fixed point 
(and ignoring higher nonlinearities) results in the nor- 
mal form of the Hopf bifurcation [l7|. The correspond- 
ing SPDE, upon ignoring noise, may be cast into the 
form of a CGLE, where a complex variable z encodes the 
densities' deviations from the internal fixed point, see 
Refs. 



we investigate the system's behavior in the vicinity of 



dtz = DAz + (ci — iuj)z — 02(1 — ics] 



(7) 
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FIG. 2: The dependence of the wavelength A on the rate 
7. Analytical predictions from the CGLE (|7| (divided by a 
factor 1.55), red line, are compared to numerical results (□). 
For 7 < 7£; 2 spirals are convectively instable. When 
approaching the bifurcation point, 7^0, computation of the 
correlation length (0) shows that the spatial structures are no 
longer determined by the (diverging) wavelength, but reach a 
constant size /c, see text. 
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FIG. 3: Regimes of stable, unstable, and neutrally stable 
biodiversity. They are revealed by computing the extinction 
probability for 7 = and p = 1, depending on D, after a wait- 
ing time t = A/", for increasing system size N. We show curves 
for = 20 X 20 (□), AT = 25 x 25 (o), AT = 40 x 40 (a), and 
A^ = 100 X 100 (0). The transitions occur at D^^^ ^ 3.5x 10~^ 



and D, 



(2) 



; 0.2. 



The coefficients cj, ci, C2, and C3 are rational functions of 
the parameter 7 and given by a; = -^(7+2), ci = ^, C2 = 

73(77'+277+27)+207^ . _ i 9^3(7+2) [237'+63(7+l)] 



73(772+277+27)+2072 



^ '90(772+277+27) ' ' ^3 

As main characteristics, near the bifurcation point, 
i.e. 7 <C 1, they asymptotically behave as const., 
ci ~ C2 ~ 7, and cs ~ I/7. The vanishing of ci and C2 at 
the bifurcation point reflects the neutral stability of the 
fixed point and the surrounding closed orbits. 

A CGLE is accurate in the vicinity of a supercritical 
Hopf-bifurcation. In our case, Eq. ([7j) is only approxi- 
mate: the RE dl]) do not exhibit a Hopf bifurcation, but 
a degenerate bifurcation where heteroclinic orbits turn 
into a family of nestled, neutrally stable orbits. In the 
following, we show that the CGLE ([7j) nevertheless pro- 
vides a reliable description of the system. However, ig- 
noring higher nonlinearities induces certain quantitative 
deviations from numerical findings, as discussed below. 

Rotating spiral waves constitute a generic solu- 
tion to the CGLE, and the spreading velocities, fre- 
quencies, and wavelengths can be calculated analyti- 
cally [2l|. As an example, the wavelength follows as 

A = 27rc3 ^JD jcx [\ — y^l + . Comparing these an- 

alytic values to numerical ones, we have found that the 
former exceed the latter by a factor of 1.55, which we at- 
tribute to higher (ignored) nonlinearities. Reducing the 
analytical result by this constant factor yields an excel- 
lent agreement, see Fig. [2l 

The CGLE ([7j) predicts an Eckhaus instability. 
Namely, the spirals are only stable against longitudi- 
nal long-wave perturbations if the Eckhaus criterion 
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l+c^ 



> is fullfilled, with the rescaled wavevector 



= 27i^D/ci/\ see Ref. In our case, this condi- 



tion translates into 7 > 7£;, with the value 7|."^^- ^ 1.43. 
Above 7£;, the spirals are absolutely stable, while for val- 
ues of 7 below 7^;, they exhibit convective instabilities: 
a localized perturbation grows but travels away. The 
instabilities result in the blurring seen in Fig. [T] and orig- 
inate in the only weak instability of the RE's internal 
fixed point. Numerically, the value of 7^; may be de- 
termined by analyzing the influence of perturbations, as 
we describe in the Supplementary Material [22[. A value 
7^; ~ 2 is found, which exceeds the analytical one at a 
factor 1.4. As for the wavelength, we attribute this devia- 
tion to the fact that the CGLE ([7j) is only a (third-order) 
approximation to the full nonlinear terms appearing in 
the SPDE (p. 

Upon approaching the bifurcation point, i.e. when 
7^0, the spirals' wavelength, as predicted by the 
CGLE ([7j), diverges: To leading order in I/7, we cal- 
culate A^^y"^P*- = 47rV'2L)/7, such that A^^^^^p*- ^ 00 
when 7^0. However, in this limit, the blurring due 
to the convective instability dominates over the instabil- 
ity of the RE (which indeed vanishes at the bifurcation 
point), such that the spiral waves are no lon ger relevant. 
Computation of spatial correlation functions [25[ and the 
resulting correlation length /corr (where spatial correla- 
tions have decayed to 1/e of their maximal value) shows 
that, for 7 < 0.001, /corr is no longer proportional to 
the (diverging) wavelength, but appears to approach a 
constant value IcVD, see Fig. [2l The spatial structures 
emerging at the bifurcation point are extremely weak, c.f. 
Fig.[T](e), and should be caused by fluctuations alone, as 
the RE do not yield an instability there. Determining 
the value Ic requires understanding of these fluctuation- 
driven patterns, and may be the subject of future studies. 

For high mobilities, the system becomes effectively 
well-mixed, and solely the stability of the RE's internal 
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fixed point determines wiietiier species diversity will be 
maintained or not |23|]. In contrast, when mobility lies 



below a threshold value, spatial patterns can form and 
help to enable stable species diversity [5[. In the follow- 
ing, we show that, at the bifurcation, patterns can have 
even different, highly non-trivial impact on biodiversity. 

Fluctuations unavoidably lead to ultimate extinction 
of species [5]. However, transient coexistence can last 
very long. For this reason, we have proposed a scheme to 
differentiate stable from unstable diversity which is based 
on time-scales. In brief, neutral stability leads to a mean 
extinction time T which is proportional to the system size 
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N 



We therefore consider the asymptotic limit 
oo, and define a situation where T/7V ^ oo as sta- 
ble diversity (extinction takes very long), while T/N 
coresponds to unstable diversity (extinction is fast). 

We have applied this concept to determine the influ- 
ence of mobility on diversity's stability at the bifurca- 
tion point 7 = 0. From stochastic simulations of the 
reactions ([T]) with nearest-neighbor exchanges on lattices 
of increasing size N ^ we have computed the extinction 
probability Pext that, starting at a random initial state 
with equal densities of A, 5, and C, two species have 
gone extinct after a waiting time t = N. Results are 
shown in Fig. [3l where three different mobility regimes 
emerge. For low mobilities (around D = 0.001), Pext ap- 
proaches for increasing implying T/N ^ oo. In this 
regime, coexistence is therefore stable, as has previously 
been observed for vanishing mobility, see fH] and Refs. 
therein. For intermediate mobilities around D — 0.007, 
the opposite behavior emerges: Pext tends to 1, such that 
T/N and coexistence is unstable. High mobilities, 
around D = 0.5, yield an asymptotic value of Pext of 
about 0.53, independent of N^ and thus neutral stabil- 
ity [Hf (the precise value of Pext depends on the choice 
of t ^ A^). Critical mobility values separate these three 
regimes. Namely, with increasing N, the transition from 
the stable to the unstable regime becomes sharp at a 
value P^c^^ ^ 3.5 x 10~^. The unstable regime verges on 
the neutrally stable one around Dc ~ 0.2, where Pext 
reaches 0.53. Our data do not reveal whether a sharp 
transition or a crossover results when N ^ oo. 

In summary, we have quantitatively analyzed the emer- 
gence of blurred (convectively instable) spiral waves near 
a bifurcation point of cyclic dynamics. This bifurcation 
point is characterized by neutrally stable, cycling orbits 
predicted by the corresponding RE (|4]). There, spa- 
tial structures are predominantly determined by noise. 
These patterns have ambiguous impact on maintaining 
the coexistence of the interacting species. Low individu- 
als' mobility, corresponding to small patterns, promotes 
diversity, as has already been observed in previous stud- 
ies [El, [gI, [ll| . However, medium values of mobility, induc- 
ing relatively large patterns, lead to rapid species extinc- 
tion. Only for high mobilities, spatial patterns have no 



influence, and, being at the bifurcation point, neutrally 
stable coexistence emerges. Further investigations of the 
destabilizing influence of spatial patterns at medium mo- 
bilities are required for a general understanding of the 
effects of spatial degrees of freedom on the coexistence of 
mobile individuals. Also, such studies will shed further 
light on the role of noise that becomes a major player in 
pattern formation at bifurcation points. 
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Supplementary Material 



Numerical determination of the onset of convective instability 



In this Supplementary Material, we describe further how the convective instability takes shape in the system's 
reactive steady state, and how one can therefrom infer numerically the value 7^;, i.e. the boundary between stable 
and convectively unstable spiral waves. 

We have already seen that convective instability manifests in a blurring of the spirals. However, at the onset of 
convective instability, the blurring is extremely weak, and cannot be detected from analyzing snapshots. Moreover, 
further inside the region of convective instability, where the blurring is stronger, it can hardly be quantified using, 
e.g., correlation functions. For these reasons, the blurring seems not to be a good candidate for a quantification of 
the onset of convective instability. 

A better suited measure is obtained by investigating the effects of convective instability on large spiral waves. To 
obtain the latter, we solve the SPDE ([5]) with very low noise, but starting from an initial state which is spatially 
inhomogeneous. The latter inhomogeneities serve as "seeds" for the developing spiral waves and determine the position 
and numbers of their vortices [13] • In our numerical solutions, we have started from initial densities a(r, t = 0) = 
a* + 1/100 X cos(27rrir2), 6(r, t = 0) = 6*, c(r, t = 0) = c*. Hereby, a*,b*,c* denote the densities at the reactive fixed 
point of the RE (gj). 



FIG. 4: Snapshots of the reactive steady state for different values of the rate 7, and a diffusivity oi D = 10~^ x 7. We have 
solved the deterministic PDF, meaning Fqs. (|5]) in the absence of noise, with initial spatial inhomogeneities: a(r,t = 0) = 
a* + 1/100 cos(27rrir2), b{r, t — 0) = b* , c(r, t — 0) — c* . From left to right, the value of 7 decreases, approaching the Fckhaus 
instability. At the same time, the spiral density rises. 

Deep in the regime of stable spiral waves, i.e. for 7 ^ 7£;, few large spirals emerge from these initial inhomogeneities, 
see Fig. |4] (a). For decreasing values of 7, approaching the convective instability, more and more vortices appear, 
see Fig. |4] (b),(c). Apparently, upon approaching the onset of convective instability, the large spiral waves become 
increasingly unstable against perturbations and noise, such that a larger number of smaller spirals appears. The density 
of spiral vortices increases, upon saturating at a constant value of about 0.3 for 7 < 2. Interestingly, the same spiral 
density also emerges when the SPDE ([5]) are solved starting without initial inhomogeneities, such that noise perturbs 
the system initially and leads to the emergence of entangled spiral waves. Consequently, at the onset of the convective 
instability (as well as in the regime of convective instability), upon starting from initial spatial inhomogeneities or 
employing noise as an internal source of inhomogeneities, the system reaches statistically equivalent states. This is 
not true within the regime of stable spirals: There, initially imposed spatial inhomogeneities can lead to very different 
vortex densities than emerge when fluctuations determine the formation of entangled spirals. 

Different vortex densities manifest in differences in the spatial correlation functions. The equal-time autocorrelation 
^aa(^) of species A at distance r shows damped oscillations, and may be approximated by 



with constants a, 6, c. Hereby, the Bessel function Jq describes oscillations, at a length proportional to 6, which relate 
to the wavelength of the spirals. The exponential damping exp(— r/6) stems from the entangled spiral waves. The 
typical length 6, at which the oscillations decay, reveals information about the vortex density: b is proportional to the 
mean distance between spiral vortices, and therefore proportional to the inverse square root of the vortex density. 




gaa{r) ^ aexp(-r/6)Jo(r/c) , 
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FIG. 5: The decay length b of the oscillations of the correlation functions versus the rate 7. We show results from solving the 
SPDE JH both starting from spatial inhomogeneities (red squares) as well as from initially homogeneous states (black circles). 
In the regime of stable spirals, the decay length is larger for the solution from initial inhomogeneities, corresponding to larger 
spirals. Approaching the Eckhaus instability, the decay length approaches the values of the initially homogeneous system. As 
a remark, similar to the wavelength, h is asymptotically (for 7 ^ 7^; as well as for 7 ^ 7^;) proportional to 7~^^^. The data 
have been obtained by numerical solution of the SPDE (|5]) with diffusivity D — 10~^ x 7. 



From numerical simulations, obtained from initial inhomogeneities, we have determined the value of h for different 
rates 7. As expected, for large values of 7, 7 ^ 7^;, i.e. in the regime of stable spirals, h is large, corresponding to 
a small vortex density arising from the small density of initially imposed spatial inhomogeneities. Decreasing 7, the 
value of h decreases, and for 7 < 7^; 2 reaches the same value as arises when stochastic effects induce spiral waves. 



